# Sensitivity Analysis 
m.out.1 <- mediate(m.comp.1, y.comp.1, treat="rel_banperm", mediator="liberal", conf.level=.95)
m.out.2 <- mediate(m.comp.2, y.comp.2, treat="rel_banperm", mediator="liberal", conf.level=.95)
m.out.3 <- mediate(m.comp.3, y.comp.3, treat="rel_banperm", mediator="liberal", conf.level=.95)
m.out.4 <- mediate(m.comp.4, y.comp.4, treat="rel_banperm", mediator="liberal", conf.level=.95)
m.out.5 <- mediate(m.comp.5, y.comp.5, treat="rel_banperm", mediator="liberal", conf.level=.95)

sen.out.1 <- medsens(m.out.1)
sen.out.2 <- medsens(m.out.2)
sen.out.3 <- medsens(m.out.3)
sen.out.4 <- medsens(m.out.4)
sen.out.5 <- medsens(m.out.5)

c.1 <- sen.out.1$d0
c.2 <- sen.out.2$d0
c.3 <- sen.out.3$d0
c.4 <- sen.out.4$d0
c.5 <- sen.out.5$d0

s.1 <- (sen.out.1$upper.d0 - sen.out.1$d0)/1.96  
s.2 <- (sen.out.2$upper.d0 - sen.out.2$d0)/1.96  
s.3 <- (sen.out.3$upper.d0 - sen.out.3$d0)/1.96  
s.4 <- (sen.out.4$upper.d0 - sen.out.4$d0)/1.96  
s.5 <- (sen.out.5$upper.d0 - sen.out.5$d0)/1.96  

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res <- round(combined.res, 2)

pdf("figures/figure_S3.pdf", width=7, height=5)
par(mfrow=c(1,2), mar=c(4,4,2,1), mgp=c(2, .5, 0), fg="black", tck=-.02, las=1)
plot(sen.out.1$rho, combined.res[,1], type="l", axes=F, xlab=expression(rho), ylab="ACME \n (Restrictive Policy Condition)",  ylim=c(-3.5, 2.5))
grid(lty=1, lwd=.5, col="lightgrey")
polygon(c(sen.out.1$rho, rev(sen.out.1$rho)), c(combined.res[,1]+1.96*combined.res[,2], rev(combined.res[,1]-1.96*combined.res[,2])), col=rgb(0, 0, 0, 80, max=255), border=F)
abline(h=0, lty=2)
abline(v=0, lty=2)
abline(h=m.out.1$d0, col="maroon3")
axis(1, col="white")
axis(2, col="white")

c.1 <- sen.out.1$d1
c.2 <- sen.out.2$d1
c.3 <- sen.out.3$d1
c.4 <- sen.out.4$d1
c.5 <- sen.out.5$d1

s.1 <- (sen.out.1$upper.d1 - sen.out.1$d1)/1.96  
s.2 <- (sen.out.2$upper.d1 - sen.out.2$d1)/1.96  
s.3 <- (sen.out.3$upper.d1 - sen.out.3$d1)/1.96  
s.4 <- (sen.out.4$upper.d1 - sen.out.4$d1)/1.96  
s.5 <- (sen.out.5$upper.d1 - sen.out.5$d1)/1.96  

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res <- round(combined.res, 2)

plot(sen.out.1$rho, combined.res[,1], type="l", axes=F, xlab=expression(rho), ylab="ACME \n (Liberal Policy Condition)", ylim=c(-3.5, 2.5))
grid(lty=1, lwd=.5, col="lightgrey")
polygon(c(sen.out.1$rho, rev(sen.out.1$rho)), c(combined.res[,1]+1.96*combined.res[,2], rev(combined.res[,1]-1.96*combined.res[,2])), col=rgb(0, 0, 0, 80, max=255), border=F)
abline(h=0, lty=2)
abline(v=0, lty=2)
abline(h=m.out.1$d1, col="maroon3")
axis(1, col="white")
axis(2, col="white")
dev.off()

pdf("figures/figure_S4.pdf", width=7, height=5)
par(mfrow=c(1,2), mar=c(4,4,2,1), mgp=c(2, .5, 0), fg="black", tck=-.02, cex.axis=.8, cex.lab=.8, ask=F, las=1)
plot(sen.out.1, sens.par = "R2", r.type = "total", sign.prod = "positive", ylab="Proportion of variance in Y \n explained  by an unobserved confounder", xlab="Proportion of variance in M \n explained  by an unobserved confounder", main="", xlim=c(0, .8), ylim=c(0, .6), axes=F, 
     ask=F)
dev.off()
